home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / cheb / eval.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.0 KB  |  178 lines
  1. /* cheb/eval.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_chebyshev.h>
  25.  
  26. /* For efficiency there are separate implementations of each of these
  27.    functions */
  28.  
  29. double
  30. gsl_cheb_eval (const gsl_cheb_series * cs, const double x)
  31. {
  32.   size_t i;
  33.   double d1 = 0.0;
  34.   double d2 = 0.0;
  35.  
  36.   double y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
  37.   double y2 = 2.0 * y;
  38.  
  39.   for (i = cs->order; i >= 1; i--)
  40.     {
  41.       double temp = d1;
  42.       d1 = y2 * d1 - d2 + cs->c[i];
  43.       d2 = temp;
  44.     }
  45.  
  46.   return y * d1 - d2 + 0.5 * cs->c[0];
  47. }
  48.  
  49. double
  50. gsl_cheb_eval_n (const gsl_cheb_series * cs, const size_t n, const double x)
  51. {
  52.   size_t i;
  53.   double d1 = 0.0;
  54.   double d2 = 0.0;
  55.  
  56.   size_t eval_order = GSL_MIN (n, cs->order);
  57.  
  58.   double y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
  59.   double y2 = 2.0 * y;
  60.  
  61.   for (i = eval_order; i >= 1; i--)
  62.     {
  63.       double temp = d1;
  64.       d1 = y2 * d1 - d2 + cs->c[i];
  65.       d2 = temp;
  66.     }
  67.  
  68.   return y * d1 - d2 + 0.5 * cs->c[0];
  69. }
  70.  
  71.  
  72. int
  73. gsl_cheb_eval_err (const gsl_cheb_series * cs, const double x,
  74.            double *result, double *abserr)
  75. {
  76.   size_t i;
  77.   double d1 = 0.0;
  78.   double d2 = 0.0;
  79.  
  80.   double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  81.   double y2 = 2.0 * y;
  82.  
  83.   double absc = 0.0;
  84.  
  85.   for (i = cs->order; i >= 1; i--)
  86.     {
  87.       double temp = d1;
  88.       d1 = y2 * d1 - d2 + cs->c[i];
  89.       d2 = temp;
  90.     }
  91.  
  92.   *result = y * d1 - d2 + 0.5 * cs->c[0];
  93.  
  94.   /* Estimate cumulative numerical error */
  95.  
  96.   for (i = 0; i <= cs->order; i++)
  97.     {
  98.       absc += fabs(cs->c[i]);
  99.     }
  100.  
  101.   /* Combine truncation error and numerical error */
  102.  
  103.   *abserr = fabs (cs->c[cs->order]) + absc * GSL_DBL_EPSILON;
  104.  
  105.   return GSL_SUCCESS;
  106. }
  107.  
  108. int
  109. gsl_cheb_eval_n_err (const gsl_cheb_series * cs,
  110.              const size_t n, const double x,
  111.              double *result, double *abserr)
  112. {
  113.   size_t i;
  114.   double d1 = 0.0;
  115.   double d2 = 0.0;
  116.  
  117.   double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  118.   double y2 = 2.0 * y;
  119.  
  120.   double absc = 0.0;
  121.  
  122.   size_t eval_order = GSL_MIN (n, cs->order);
  123.  
  124.   for (i = eval_order; i >= 1; i--)
  125.     {
  126.       double temp = d1;
  127.       d1 = y2 * d1 - d2 + cs->c[i];
  128.       d2 = temp;
  129.     }
  130.  
  131.   *result = y * d1 - d2 + 0.5 * cs->c[0];
  132.  
  133.   /* Estimate cumulative numerical error */
  134.  
  135.   for (i = 0; i <= cs->order; i++)
  136.     {
  137.       absc += fabs(cs->c[i]);
  138.     }
  139.  
  140.   /* Combine truncation error and numerical error */
  141.  
  142.   *abserr = fabs (cs->c[cs->order]) + absc * GSL_DBL_EPSILON;
  143.  
  144.   return GSL_SUCCESS;
  145. }
  146.  
  147. int
  148. gsl_cheb_eval_mode_e (const gsl_cheb_series * cs,
  149.               const double x, gsl_mode_t mode,
  150.               double *result, double *abserr)
  151. {
  152.   size_t i;
  153.   double d1 = 0.0;
  154.   double d2 = 0.0;
  155.  
  156.   double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  157.   double y2 = 2.0 * y;
  158.  
  159.   size_t eval_order;
  160.  
  161.   if (GSL_MODE_PREC (mode) == GSL_PREC_DOUBLE)
  162.     eval_order = cs->order;
  163.   else
  164.     eval_order = cs->order_sp;
  165.  
  166.   for (i = eval_order; i >= 1; i--)
  167.     {
  168.       double temp = d1;
  169.       d1 = y2 * d1 - d2 + cs->c[i];
  170.       d2 = temp;
  171.     }
  172.  
  173.   *result = y * d1 - d2 + 0.5 * cs->c[0];
  174.   *abserr = fabs (*result) * GSL_DBL_EPSILON + fabs (cs->c[eval_order]);
  175.  
  176.   return GSL_SUCCESS;
  177. }
  178.